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piji Abstract 

Ph ■ Biological functions are generated as a result of developmental dynam- 

O ! ics that form phenotypes governed by genotypes. The dynamical system for 

o I development is shaped through genetic evolution following natural selection 

qh' based on the fitness of the phenotype. Here we study how this dynamical 

'""'■ system is robust to noise during development and to genetic change by muta- 

,-H ■ tion. We adopt a simplified transcription regulation network model to govern 

^ ! gene expression, which gives a fitness function. Through simulations of the 

\Q I network that undergoes mutation and selection, we show that a certain level 

cn ■ of noise in gene expression is required for the network to acquire both types 

of robustness. The results reveal how the noise that cells encounter during 

^ . development shapes any network's robustness, not only to noise but also to 

OO ! mutations. We also establish a relationship between developmental and mu- 

.. [ tational robustness through phenotypic variances caused by genetic variation 

.^ ■ and epigenetic noise. A universal relationship between the two variances is 

^ ' derived, akin to the fiuctuation-dissipation relationship known in physics. 

^ '. Lead paragraph 

Biological function is generally robust to noise in dynamical sys- 
tems and to mutation to structure of the system. How are these 
two types of robustness related? Through numerical simulations 
of a minimal model that captures the essence of a gene expres- 
sion dynamics that undergoes a mutation and selection process, 
we demonstrate that the evolved networks acquire robustness to 
mutation only when gene expression is sufficiently noisy - a physical 
constraint highlighted by a series of experimental studies in recent 



years. We derive a quantitative condition for evolution of robust- 
ness which ultimately links robustness to mutation in the evolu- 
tionary time-scale and robustness to noise in development in the 
reproductive time-scale. Our model simulations predict correla- 
tion between the two robustness, which lead to a universal propor- 
tionality relationship between variances of genetic and epigenetic 
phenotype fluctuations, which remind us of fluctuation-dissipation 
relationship in statistical physics. 

1 Introduction 

Function in a biological system is generated by dynamical systems in gen- 
eral. A specific shape of a protein that gives rise to some function is formed 
by folding dynamics. In a cell, gene expression pattern changes in time, to 
shape some chemical composition, from which a given function of the cell 
emerges. In multicellular organisms, developmental dynamics lead to pat- 
terns of differentiated cells. Phenotype is shaped in a broad sense from such 
'developmental dynamical systems', according to which the fitness of an in- 
dividual is determined. 

In a biological system, such developmental dynamics are shaped through 
Darwinian evolution. Offspring are produced depending on the fitness of 
the parents, with some variations that introduce slight modifications in pa- 
rameters or networks in developmental dynamical systems. Among modified 
dynamical systems, phenotypes with higher function are selected for the next 
generation. Hence, evolutionary shaping of function is a result of variation 
and selection of developmental dynamical systems. 

Combining both the terminology of biology and dynamical systems, evo- 
lutionary process is summarized as follows, (i) The genotype gives a set of 
equations, i.e., terms and parameters in developmental dynamical systems, 
(ii) From an initial condition, the developmental dynamical system leads to a 
set of state variables for a phenotype. (iii) Fitness is a function of these state 
variables. The offspring number increases with fitness, (iv) In reproduction, 
there are mutations (or other sources for variations) giving rise to variation 
in dynamical systems. At each generation, there is redistribution of dynam- 
ical systems. If evolution serves to increase fitness, dynamical systems with 
higher function are shaped through this selection-mutation process. 



In discussing this evolutionary shaping of dynamical systems with higher 
function, an important issue is robustness. Robustness is defined as the 
ability to function against changes in the parameters of a given system [1, 
2, 3, 4, 5, 6]. In any biological system, these changes have two distinct 
origins: genetic and epigenetic. The former concerns structural robustness 
of the phenotype, i.e., rigidity of phenotype against variation in dynamical 
systems, introduced by genetic changes produced by mutations. On the 
other hand, the latter concerns robustness against the stochasticity that can 
arise in a given dynamical system, which includes fluctuation in initial states 
and stochasticity occurring during developmental dynamics or in the external 
environment. For example, stochasticity in gene expression has recently been 
studied extensively both experimentally and theoretically [7, 8, 9, 10, 11]. 
Indeed, the existence of such stochasticity is natural, because the number 
of molecules in a cell is generally limited. Accordingly, phenotypes of cells 
differ, even among those sharing the same genotype [12]. An important recent 
recognition is that phenotypic noise is indeed significant, as is highlighted by 
the log- normal distribution of protein abundances in bacteria[9]. 

The two types of robustness — epigenetic and genetic — are concerned with 
fluctuations in phenotype caused by noise in developmental processes and by 
variations in dynamical systems caused by genetic changes. In terms of dy- 
namical systems, these two types of robustness are the stability of a state (an 
attractor) to external noise and the structural stability of the state against 
changes in the underlying equations, respectively. Then, how are these two 
forms of robustness related in evolution? Recently we have found a possible 
relationship between the two, from an experiment on adaptive evolution in 
bacteria[13] and numerical evolution of reaction networks and transcription 
regulation networks [5, 14]. The two types of robustness, measured by the 
(inverse of) phenotypic variances, increase in correlation through evolution, 
whereas evolutionary speed is greater as the magnitude of the phenotypic 
fluctuation due to noise is increased. Hence, the relevance of noise to evolu- 
tion is suggested. 

Despite recent quantitative observations on phenotypic fluctuation, noise 
is often thought to be an obstacle in tuning a system to achieve and maintain 
a state with higher functions, because the phenotype may be deviated from 
an optimal state achieving a higher function through such noise. Indeed, the 
question most often asked is how a given biological function can be main- 
tained in spite of the existence of phenotypic noise[ll, 15]. Although the 



positive roles of stochasticity in gene expression to cell differentiation [16] 
and adaptation have been discussed [17, 18], its role in evolution has not 
been explored fully. As a relatively large amount of phenotypic noise has 
been preserved through evolution, it is important to investigate any positive 
roles of such noise for the evolution of biological functions. 

The present paper is organized as follows. Following the Ref. [5], we intro- 
duce a simple transcription regulation network model in §2. In this model, 
fitness is determined by the gene expression pattern generated by this net- 
work. Using a genetic algorithm to change the network so that the fitness 
is increased, we show that evolution of the robustness of fitness to muta- 
tion to a given network is possible only under a certain level of noise in 
the gene expression dynamics. Dependence of the fitness distribution on the 
noise level is analyzed in §3. Variances of fitness due to developmental noise 
and to genetic variation are computed, which give indices for developmen- 
tal and mutational robustness, respectively. The proportional relationship 
between the two variances is given in §4, whereas in §5, the two variances 
are computed for the expression of all genes, which demonstrate common 
proportion relationship, suggesting the existence of a universal relationship 
akin to fluctuation-dissipation theorem in statistical physics. Summary and 
discussion are given in §6. 

2 Model 

To study evolutionary shaping of dynamical systems, we consider a simple 
model for 'development'. This consists of a complex dynamic process to 
reach a target phenotype under a condition of noise that might alter the 
final phenotypic state. We do not choose a biologically realistic model that 
describes a specific developmental process, but instead take a model as simple 
as possible to satisfy the minimal requirement for our study. Here we take 
a simplified model for gene expression dynamics, governed by transcription 
regulatory networks. Expression of a gene activates or inhibits expression of 
other genes under noise. These interactions between genes are determined by 
the transcription regulatory network (TRN). The expression profile changes 
in time, and eventually reaches a stationary pattern. This gene expression 
pattern determines fitness. 

To be specific, a typical switch-like dynamics with a sigmoid input-output 



behavior [19, 20, 21] was adopted, although several simulations in the form 
of biological networks will give essentially the same result. In this simplified 
model, the dynamics of a given gene expression level Xi is described by: 
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dxi/dt = tanh[^=^=== V JijXj] - Xi + (Jr]i(t), (1) 

VM-k ^^^ 

where Jjj = —1,1, 0, and rii{t) is Gaussian white noise given by < rii(t)r]j{t') >= 
6ij6{t — t'). M is the total number of genes, and k is the number of out- 
put genes that are responsible for fitness to be determined. The value of a 
represents noise strength that determines stochasticity in gene expression. 
By following a sigmoid function tank, Xi has a tendency to approach either 
1 or -1, which is regarded as 'on' or 'off' in terms of gene expression. The 
initial condition is given by (-1,-1,. ..,-1); i.e., all genes are off unless noted 
otherwise. 

Depending on the expression pattern, the cell state changes. Hence, the 
function of the cell and accordingly its fitness changes. Here we assume that 
fitness is determined by setting a target gene expression pattern. In the 
present paper we adopt a simple target so that the gene expression levels 
(xj) for the output genes i = 1,2,- ■ ■ ,k < M reach 'on' states, i.e., Xi > 0. 
The fitness F is at its maximum if all k genes are 'on' after a transient time 
span Tini, and at its minimum if all are off. F is set at 0, if all the target 
genes are on, and is decreased by 1 if one of the k genes is 'off'. Accordingly, 
the fitness function is defined by 

F = ^ Y. C (Stgnix,) - l)dt, (2) 

^y-'- f -'-ini) j—\JTini 

where Sign{x) = 1 for x > and —1 otherwise. The initial time Tj„j can 
be considered as the time required for developmental dynamics. The fitness 
is computed only after time Tj„j, which is sufficiently large for a given gene 
expression's dynamics to fall on an attractor. In the simulations presented 
here we adopt Tj„j = 50 and the results to be discussed are not altered even 
if it is increased. Here, in considering also the possibility of an oscillatory 
attractor, we adopt the temporal average for the fitness, but indeed, it is not 
necessary as attractors after evolution are mostly fixed points. Tj is set at 



Tini + 50 here, but even if it is defined just by a snapshot value at t = Tj„j, 
the results are not essentially altered. 

Selection is applied after the introduction of mutation at each generation 
in the TRN. Among the mutated networks, we select those with higher fitness 
values. Because the network is governed by Jij which determines the 'rule' of 
the dynamics, it is natural to treat Jij as a measure of genotype. Individuals 
with different genotype have a different set of Jij. 

As model (1) contains a noise term, the fitness can fluctuate at each run, 
which leads to a distribution in F, even among individuals sharing the same 
network. For each individual network, we compute the average fitness F 
over a given number of runs. At each generation there are A^ individuals 
with different sets of Jij. Then we select the top Ns{< N) networks that 
have higher fitness values. 

At each generation there are A^ individuals. We compute the average 
fitness F for each network by carrying out L runs for each. Then, A^^ = N/4 
networks with higher values of F are selected for the next generation, from 
which Jij is 'mutated', i.e., Jij for a certain pair i,j selected randomly with 
a certain fraction is changed among ±1, 0. (To be specific, only the changes 
f ^^ ^^ — 1 are allowed). The fraction of path change is given by the 
mutation rate /i. Unless otherwise mentioned, only a single path, i.e., a single 
pair of i,j is changed so that the mutation rate n = 1/M{M — k). Here we 
make N/Ng mutants from each of the top Ns networks, to keep A^ networks 
again for the next generation. Following mutation, the A^ individuals at 
each generation have slightly different network elements, Jij, so that the 
values of F differ. From this population of networks we repeat the processes 
of developmental dynamics, their mutation, and selection of networks with 
higher fitness values. 

Unless otherwise mentioned, we chose N = L = 200, while the conclusion 
to be shown does not change as long as these are sufficiently large. (We have 
also carried out the selection process by F instead of F, but the conclusion 
is not altered if A^ is chosen to be sufficiently large.) Throughout the paper, 
we use /5 = 50. For most numerical results we use M = 64 and k = 8, unless 
otherwise mentioned. Initially we chose Jij randomly with equal probability 
for ±1, 0. However, the results to be discussed are not changed qualitatively, 
even if the fraction of is increased. As shown in eq.(l), we did not include 
a connection from the output genes i < k, so that Jij with j < k is fixed at 
and is not mutated. However, the results to be discussed are not altered, 

6 



even if such connections are included in eq.(l). 

3 Decrease in the average fitness with the de- 
crease of noise 

Let us first see how the evolutionary process changes as a function of the 
noise strength a. Within a hundred generations, the top fitness among the 
network population approaches the highest value 0, if noise level is not too 
high {a < a'^ ^ .12 for M = 64) (see Fig. 1). However, this is not the 
case for the average fitness among a population with slightly different genes. 
For low noise level, average fitness stays lower than the fittest value. In the 
middle range of noise level, the average value approaches the fittest value. 
This difference against the noise level is clearer in the temporal evolution 
of the lowest fitness among the existing networks. As shown in Fig. 1, the 
value stays rather low for a small noise case, whereas for evolution under 
a high noise value, it goes up with the top fitness value. In other words, 
the distribution of the fitness P{F) over existing networks has three distinct 
behaviors, depending on the noise strength a. 

(i) For small o" (< (Jc ~ 0.04 for M = 64 and the mutation rate with 1 
path per generation), the distribution is broad. The top reaches the fittest, 
although there remain individuals with very low fitness values F, even after 
many generations of evolution. 

(ii) For the middle range a {cc < cr < a^), the distribution is sharp and 
concentrated at the fittest value. Even those individuals with the lowest 
fitness approach F = . We call this range the 'robust evolution region'. 

(iii) For the much larger noise value {a > a'J, the top does not reach the 
fittest values, while the distribution is sharp, and concentrated around the 
top value. Here the noise is so large that Xi changes its sign at a certain 
fraction even after reached the target, so that the top fittest is not achieved. 

The change in the distribution between (i) and (ii) is demonstrated in 
Fig. 2. There is a threshold noise o"c, below which the distribution P{F) 
is broadened. As a result, the average fitness over all individuals, < F >= 
J FP{F)dF is low. < F > and the lowest fitness values over individuals Fmin, 
after a sufficiently large number of generations, are plotted against a in Figure 
3. The abrupt decrease in fitness suggests the existence of threshold noise 
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Figure 1: Evolutionary time course of the fitness F. The highest, average, 
and lowest values of the fitness F an&ng all individuals that have different 
genotypes (i. e., networks J^) at each generation are plotted, (a) a = 0.01, 
(b) a = 0.1 (c)or = 0.2. 
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Figure 2: Fitness distribution P{F). From the top fitness network evolved 
after 200 generations, we generated 1000 networks by clianging a single el- 
ement in the Ji,jS matrix, and computed the average fitness F for each, to 
obtain the fitness distribution. Inset is the magnification for —0.2 < F < 0. 
The histogram is computed with a bin size 0,01, whereas for the inset, we 
adopt the bin size 0.001. For high a (red, with a = 0.1), the distribution is 
concentrated at F = 0, whereas for low a (green, with a = 0.01), the distri- 
bution is extended to large negative values, even after many generations. 



level (Tc, below which low-fitness mutants always remain in the distribution. 

This transition on the noise level is rather sharp, as long as the number of 
genes M is large. As M becomes smaller, the plateau in the highest fitness 
value against a is decreased, and the transition loses sharpness, as shown 
in Fig. 4 for M = 16, with the mutation rate per single path change per 
generation. With an increase in M, the region achieving the highest fitness 
with sharp distribution increases, as well as the increase in the sharpness of 
transition. 

If there are more genes, change in a single path has a smaller infiuence, so 
that an increase in the robust evolution region is to be expected. However, 
the increase in the robust region with M is not solely caused by the decrease 
in the effective mutation rate. In Fig. 5, we have changed the number of 
genes M, by fixing the mutation rate per path, i.e., a single path change per 
generation for M = 16, 4 paths for M = 32, and 16 for M = 64. Even in 
this fixed mutation rate, the transition is sharper, and the increase in the 
robustness in the evolution is detected. 

Of course, the increase in the mutation rate reduces the robust evolution 
region. Comparing Fig. 2(a) and Fig. 4 shows that the threshold noise level 
(Tc for the robust evolution region is lowered with the decrease in mutation 
rate, where the mutation rate is 16 times higher for the simulation for M = 64 
in Fig. 4 than that for Fig. 2(a). 

At the transition at a ~ o"c we discuss here, mutational robustness of 
the evolved network changes. The top fitness network evolved at a low noise 
level {a < Cc) does not have robustness against mutation, in contrast to 
those evolved at a higher noise level. To check this distinction statistically, 
we measured the average fitness of mutants generated from the top fitness. 
We took a network with the top fitness evolved after generations, and changed 
m paths randomly (i.e., change in a value of Jij among ±1,0). By making 
a thousand networks of such m mutations, we have measured the average 
fitness of such mutants, which is plotted against m in Fig. 6. For a < 0"^, 
the average fitness decreases linearly with the number of added mutations 
m, as < F >= —C{a)m. The coefficient C{a) decreases with the increase 
in a, and for a = .1(> (Xc), the linear decrease component vanishes, giving 
rise to a plateau around < F >= 0. In other words, the fitness landscape 
has an almost neutral region. The fitness is rather insensitive to mutation, 
demonstrating the evolution of mutational robustness. Around a ~ (Xc, C{a) 
estimated at small mutation number m is rather small, and we may expect 
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Figure 3: Average of average fitness < F > and minimal fitness over gen- 
erations, plotted against the noise strength a. < F >, the average of the 
average fitness F over all individuals is computed for 100-200 generations. 
(For a > .1, we choose the average over 200-300 generation, as the fitness 
plateau is reached later). The minimal fitness is computed from the time 
average of the least fit network present at each generation, (a) M = 64 
(b)M = 16. The mutation rate is one path over all paths. In (b), since there 
are 16 times paths of (a), the mutation rate for each path is 1/16, compared 
from (a). 
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Figure 4: Average of average fitness < F > over generations after transients, 
plotted against tlie noise strength a. < F >, the average of the average fitness 
F over all individuals is computed over 100-200 generations for a < .1, and 
200-300 generations for a > .1 to remove the transient effect. The fitness for 
M = 16, 32, and 64 are computed, by fixing the mutation rate at 1 per 16 x 
16 paths, i.e., one path in total for M = 16, four paths for M = 32, and 16 
paths for M = 64. 
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Figure 5: Decline of the average fitness plotted as a function of m, which is 
the number of mutated paths from a top-fitness network. The average fitness 
is computed from F's over 1000 mutants generated from an evolved network 
having the top fitness, by inserting m mutations. For a < a^ it decreases 
linearly as < F >= —C{a)m, as a function of the number of mutated paths 
M. Here C{a) decreases. At a = .1 > ac, C{a) = 0. 



that C{a) ^ as a ^ cXc, although careful examination of the slope at 
m —>■ limit is needed to confirm it. 

The transition with regards to the mutational robustness at ac reminds 
us of the error catastrophe introduced by Eigen and Schuster [22], since the 
population of non-fit mutants starts to accumulate after this transition. In 
their error catastrophe, the transition occurs when the advantage by fitted in- 
dividuals cannot overcome combinatorial diversity of non-fit mutants. Their 
error catastrophe occurs with the increase in the mutation rate, while our 
robustness transition occurs with the decrease in developmental noise. In our 
case, the transition is of dynamic origin, as will be discussed in the next sec- 
tion, whereas there exist (combinatorially) much more non-fit networks than 
the fittest ones. The robustness transition here could be related with (or 
regarded as an extension of) error catastrophe, but the relationship between 
the two remains to be clarified. 
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4 Dynamic origin of the robustness transi- 
tion 

Why does the system not maintain the highest fitness state under a condition 
of small phenotypic noise a < a^'? Indeed, the top fitness networks that 
evolved under low noise have dynamic behaviors distinct from those that 
evolved under high noise. First, the highest fitness network evolved at low 
a often fails to reach the target if simulated under a higher noise level. The 
expression level often exhibits a few oscillations before reaching the target 
state and noise might cause expression of the output genes to switch to 'off' 
states. In contrast, the temporal course of gene expression evolved for a > cTc 
is much smoother, and is not affected by noise. This distinction is confirmed 
by simulating gene expression dynamics by cutting off the noise term over 
a variety of initial gene expression conditions, and to check if the orbit is 
attracted to the original target. We found that, for networks evolved under 
o" > o"c, a large portion of the initial conditions is attracted into the target 
pattern, while for those evolved under a < ac, only a tiny fraction (i.e., the 
vicinity of all 'off' states) is attracted to the target (see Fig. 6). 

If the time course of a given gene expression to reach its final pattern is 
represented as a motion falling along a potential valley, our results suggest 
that the potential landscape becomes smoother and simpler through evolu- 
tion and loses ruggedness after a few hundred generations. This 'develop- 
mental' landscape is displayed schematically in Fig. 7. For networks evolved 
under a > o"c there is a large, smooth attraction to the target, whereas for the 
dynamics evolved under a < (Xc, the initial states are split into small basins, 
from each of which the gene expression patterns reach different steady states. 

Now, consider the mutation to a network. Any change in the network 
leads to slight alterations in gene expression dynamics. In smooth dynamics, 
as in Fig. 7(a), this perturbation influences the attraction to the target only 
slightly. By contrast, under the dynamics as shown in Fig. 7(b), a slight 
change easily destroys the attraction to the target attractor. For this latter 
case, the fitness of mutant networks is distributed down to lower values, 
which explains the behavior observed in Fig. 3. 

Accordingly, evolution to eliminate ruggedness in developmental potential 
is possible only for sufficient noise amplitude, whereas ruggedness remains 
for small noise values and the developmental dynamics often fail to reach the 
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Figure 6: Distribution of the fitness value when the initial condition for Xj is 
not fixed at -1, but is distributed over [—1, 1]. We chose the evolved network 
as in Fig. 5; for each network we took 10000 initial conditions, and simulated 
the dynamics (1) without noise, to measure the fitness value F after the 
system reached an attractor (as the temporal average 400 < t < 500). The 
histogram is plotted with a bin size 0.1 using a semi- log plot. 
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Figure 7: Schematic representation of the basin structure, represented as 
a process of chmbing down a potential landscape. A smooth landscape is 
evolved under high noise (above), and a rugged landscape is evolved under 
low noise (below). 
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Figure 8: Relationship between Vg and V^p. Vg is computed from P{F) at 
each generation, and Vip by computing the variances of isogenic fluctuation 
of the fitness (over L runs) for existing individuals, and averaging them. 
Plotted points are over 200 generations. For a > ac ^ .03, both decrease 
with generations. See text for the definition of Vig and Vg. 



target, either by noise in gene expression dynamics or by mutations to the 
networks. It is interesting to note that a greater set of initial conditions is 
attracted to a target pattern for networks evolved under conditions of high 
noise. Existence of such global attraction in an actual gene network has 
recently been reported for the yeast cell cycle[23]. 



5 Evolution of robustness 

According to our results in the two previous sections, a network evolved under 
high noise has two types of robustness. The target pattern is reached even if 
the developmental dynamics are perturbed by noise or by a mutation to the 
network. The fitness is rather insensitive to both of such perturbations. 

As an index for robustness, variance of the fitness (or phenotype) is 
useful[5]. There are two types of variance corresponding to the above two 
types of robustness, i.e., genetic and epigenetic robustness. Variance corre- 
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spending to genetic change, denoted by Vg, is defined as phenotype variance 
caused by the distribution of genes; 

Vg = Jp(F)(F- < F >fdF, (3) 

where we note P{F) is the distribution of average fitness over all networks 
at each generation, and < F >= P{F)FdF is the average of average fitness 
over all networks. 

As the variance decreases, the system increases its robustness to genetic 
change (mutation). On the other hand, epigenetic robustness is measured by 
the phenotypic variance of isogenic organisms Vip. Indeed, the phenotypes of 
isogenic individual organisms fluctuate, as discussed extensively [7, 8, 9, 10, 
11]. We define the isogenic phenotypic variance Vip{I) by 

Vp{I) = jp{F-I){F-FifdF, (4) 

where p(-F; /) is the fitness distribution among isogenic individuals / sharing 
the same network Jjj, and Fi = J Fp{F; I)dF is the average fitness of the 
genotype J. Vip{I) generally depends on the individual (genotype) J. As a 
measure of Vip, we adopted such I that gives the peak in the fitness distri- 
bution P{F), i.e., most typical genotype. Or, instead, we estimated it by 
the average of Vip{I) over all genotypes / existing at each generation. The 
overall Vip-Vg relationship is not altered between the two estimates. In Fig.8, 
we adopted the latter estimate, to reduce statistical error. 

Under high noise conditions, the selection process favors a developmen- 
tal process that is robust against it. This robustness to noise is then em- 
bedded into robustness against mutation. Indeed, for a > a^, both Vg 
and Vip decrease through the course of evolution (Fig. 8), while main- 
taining proportionality between the two. Such proportionality between the 
two has been discussed from an evolutionary stability analysis under a few 
assumptions[14, 24, 25]. For a > ac, the the inequality Vip > V^ is satisfied, 
while it is broken at the transition a ~ cxc, which also agrees with the theo- 
retical analysis. On the other hand, the proportionality between Vip and Vg 
is also consistent with an observation on an experiment in bacterial evolu- 
tion [13], if the Fisher's theorem[26, 27, 28] on the proportionality between 
evolution speed and Vg is apphed. 
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Figure 9: The evolutionary change oiVg{i) plotted as a function of generation. 
At each generation, Vg{i) is computed from the distribution of A^ =200 indi- 
viduals at each generation, from the average value of Sign{xi) over L =200 
runs. The black lines are the variance for target genes i = 1,2, ...,k, while the 
others are for non-target genes with i = k + 1, .., M. The variance values of 
target genes decrease at earlier generations, while about a half of non-target 
genes also exhibit decreases, a = 0.1. 



6 Consolidation of the Expression of Non- 
target Genes 

So far we have studied how a dynamical system consisting of positive and 
negative feedbacks is shaped to generate a given target output pattern. Here, 
the matrix Jij has a large number of nonzero elements, and generally the 
structure is complex. Fitness generated from such a complex network is 
robust to noise and to mutation. How is such a robust system achieved to 
fixate the expressions of target genes? 

In our model, there are more degrees of freedom (genes) Xj besides those 
for the target genes. The expression levels of the non-target genes do not 
influence the fitness, and can take either positive or negative values. Hence, 
there is no selection pressure leading them to be fixed, or robust to noise or 
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Figure 10: Relationship between Vg{i) and Vip{i). Vg{i) is computed as a 
variance of the distribution of {Sign{xi)) over 200 individuals at each gen- 
eration, and Vip{i) as that of the distribution of Sign{xi) over 200 runs, 
cr = 0.1. (a) Each point is a plot from one generation for four different genes 
from one target + and three non-target (x,*,n). Plotted over 400 genera- 
tions, where (V^(i), Vip(i)) at later generations take smaller values, (b) The 
plot of (V^(^), Vip{i)) for all genes i at the generation 10(red +), 20 (green x), 
40(blue *), 80(pink D, 160(sky blue *), and 320 (black o). 
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to mutation. However, we have found that expression levels of many (about 
half) non-target genes are fixed to positive or negative values successively 
through the course of evolution in the robust evolution region, i.e., for a > a^.. 
Even after mutation to Jjj values, the sign value of Xi does not change for 
most genes i(data not shown). For such genes, the variance of each expression 
level over mutants becomes rather small because of evolution. 

To check for a possible decrease in the variance of expression levels, we 
define Vg{i) for gene i, as the variance of Sign{xi) computed at t = T^j, (i.e., 
after the time span for development) over individuals with different genes. 



instead of the variance of fitness we measured in the last section. {Sign{xi) 
is the average among isogenic individuals, i.e., the average over L runs). In 
other words, Vg{i) is defined by replacing F by Sign{xi) in eq.(3). This 
variance Vg{i) is plotted in Fig. 9 as a function of generation. Here, Vg{i) for 
i < k, the variance of the expression level of the target genes decreases at an 
earlier generation, together with a few other genes, followed by a decrease 
in Vg{i) for about half of the others towards low levels. Indeed, after 200 
generations of evolution, the variance Vg{i) for about half of the genes is 
less than 10~^. Expression levels of many degrees of freedom other than the 
requested degrees of freedom as a target are also fixed. This consolidation 
is relevant to achieve a robust system. Indeed, for a < (Xc, such fixation of 
non-target gene expression levels is not observed. 

Similarly to Vg{i), we define Vip{i) as the isogenic variance for each Sign{xi) 
by using Sign{xi) instead of F in eq. (4), and have computed its evolution 
over generations. For genes which show the decrease in Vg{i), the variances 
Vip{i) also decrease through evolution. Furthermore, for such genes i, Vg{i) 
and Vip{i) decrease in proportion, as shown in Fig. 10(a). The proportion- 
ality relationship between Vg{i) and Vip{i) holds true not only for fitness but 
also for many genes, including non-target ones. Surprisingly, as evolution 
progresses, the plot of (Vipi^i), Vg{i)) approaches a single fine, suggesting that 
the proportion coefficient between the two approaches the same value for all 
genes, whose gene expressions are fixated to on or off by generation. In fact, 
we have plotted (V^p(i), V^(i)) for all genes in Fig. 10(b). For each given 
generation, the overlaid plot of (yip{i),Vg{i)) fits on the same line for most 
genes. In other words, the proportion relationship between the two variances 
holds not only for each gene's expression through the course of evolution, but 
also for expression levels over different genes at each generation. Over gen- 
erations, the plots shift toward smaller values {Vip{i) , Vg{i)) , approaching a 
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unique line. This 'universal' coefficient for all genes suggests the existence of 
a global potential dynamic system governing many gene expression patterns. 

7 Summary and Discussion 

We have shown here that a dynamical system that is robust both to noise and 
to structural variation is shaped through evolution under noise. In our study, 
robustness to developmental noise and to mutation are represented quanti- 
tatively in terms of the phenotypic variance of isogenic organisms Vip and 
by genetic variation Vg. The proportionality between the two through evo- 
lution is a quantitative manifestation of how developmental and mutational 
robustness can evolve in coordination. In fact, whether these two types of 
robustness emerge under natural selection has long been debated in the con- 
text of developmental dynamics and evolution theory[l, 2, 29, 30], since the 
propositions of stabilization selection by Schmalhausen[31] and canalization 
by Waddington[32, 33] more than half a century ago. Here we have demon- 
strated a correlation between developmental robustness to noise and genetic 
robustness to mutation, and have shown that the former leads to the latter. 

Isogenic phenotypic fluctuation is related to phenotypic plasticity, which 
is a degree of phenotype change in a different environment. Positive roles 
for phenotypic plasticity in evolution have been discussed elsewhere [30, 34, 
35, 36]. Because susceptibility to environmental changes and phenotypic 
fluctuation are correlated positively according to the fluctuation-response 
relationship [24], our present results on the relationship between phenotypic 
fluctuations and evolution imply a relationship between phenotypic plasticity 
and evolution akin to the genetic assimilation proposed by Waddington[32]. 

Although we have demonstrated this evolution of robustness using a net- 
work model of transcriptional regulation, we expect this behavior to be ob- 
servable generally if fitness is determined through developmental dynamics 
that are sufficiently complex so that a given developmental process, when 
deviated by noise, may fail to reach the fittest target pattern. In fact, the 
decrease in phenotypic variance as well as the proportionality law between 
Vg and Vip has been con&med in a catalytic reaction network model [14]. 

We have also found that the expression of many non-target genes is also 
fixed. The behavior of many other degrees of freedom is consolidated with 
that directly related to fitness. With this consolidation, the system's be- 
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havior is constrained to exhibit an identical phenotype against noise and 
mutation. Many degrees of freedom are highly correlated, which possibly 
results in the emergence of a simple set of global dynamics over many gene 
expressions. The existence of a universal proportionality between Vg and Vip 
with an identical proportion coefficient over many genes might be a manifes- 
tation of such global dynamics. It is interesting to note that this universal 
proportionality over components is also observed in a reaction network model 
for a reproducing cell[37]. 

During the course of evolution, the variance levels of some gene expres- 
sions decrease at an early generation, while others decrease only later. This 
difference between genes should be related to how their expressions influence 
those for the target. It will be important to study the order of fixation in 
terms of the network structure. 

Borrowing the concept from statistical physics, we previously proposed 
an evolutionary fluctuation-response relationship, where proportionality be- 
tween evolution speed (response) and isogenic fluctuations is proposed and 
experimentally verifled. As the speed of evolution is proportional to Vg ac- 
cording to Fisher's fundamental theorem of natural selection[26, 27], the rela- 
tionship then suggested the proportionality between Vip and Vg. Note that in 
physics, the proportion coefficient between response and fluctuation is repre- 
sented universally in terms of temperature, and is known as the fluctuation- 
dissipation relation[38, 39] or Einstein's relation[40]. Discovery of our uni- 
versal proportion coefficient over genes suggests the existence of a general 
theory akin to the fluctuation-dissipation theorem that could be applied to 
evolutionary biology and also to robust design of dynamical systems. 

In the present paper, we adopted a simple fltness condition favoring a flxed 
expression of target genes. Accordingly, such system that has a large basin 
for a flxed point attractor corresponding to the target pattern is evolved. 
It should be of importance to examine a case with a more complex fltness 
condition. A straightforward extension is the use of several target patterns 
depending on some inputs applied to some genes. Another extension is the 
use of fltness favoring for dynamic attractors. As the model (1), with suitable 
choice of Jij, shows periodic or chaotic attractors [41], study on the evolution 
of robust dynamic system is also of importance ( see also [42]). 

I would like to thank Chikara Furusawa, Koichi Fujimoto, Masashi Tachikawa, 
Shuji Ishihara, Tetsuya Yomo, and Satoshi Sawai for stimulating discussion. 
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